AN EFFICIENT ALGORITHM FOR ACCELERATING THE CONVERGENCE 
OF OSCILLATORY SERIES, USEFUL FOR COMPUTING THE 
POLYLOGARITHM AND HURWITZ ZETA FUNCTIONS 
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ABSTRACT. This paper sketches a technique for improving the rate of convergence of 
a general oscillatory sequence, and then applies this series acceleration algorithm to the 
polylogarithm and the Hurwitz zeta function. As such, it may be taken as an extension of 
the techniques given by Borwein's "An efficient algorithm for computing the Riemann zeta 
function"(4]|l], to more general series. The algorithm provides a rapid means of evaluating 
Li,(z) for general values of complex s and a kidney-shaped region of complex z. values 
given 

by \z 2 /(z-l)\ <4. By using the duplication formula and the inversion formula, the 
range of convergence for the polylogarithm may be extended to the entire complex z-plane, 
and so the algorithms described here allow for the evaluation of the polylogarithm for all 
complex .s and z, values. 

Alternatively, the Hurwitz zeta can be very rapidly evaluated by means of an Euler- 
Maclaurin series. The polylogarithm and the Hurwitz zeta are related, in that two evalua- 
tions of the one can be used to obtain a value of the other; thus, either algorithm can be used 
to evaluate either function. The Euler-Maclaurin series is a clear performance winner for 
the Hurwitz zeta, while the Borwein algorithm is superior for evaluating the polylogarithm 
in the kidney-shaped region. Both algorithms are superior to the simple Taylor's series or 
direct summation. 

The primary, concrete result of this paper is an algorithm allows the exploration of the 
Hurwitz zeta in the critical strip, where fast algorithms are otherwise unavailable. 
A discussion of the monodromy group of the polylogarithm is included. 



1. Introduction 

This note sketches a technique for improving the rate of convergence of a general oscil- 
latory series, and then applies this technique to the computation of the polylogarithm, and, 
in particular, to the Hurwitz zeta function. It essentially generalizes an algorithm given by 
Peter Borwein for computing the Riemann zeta function. The principle result is an algo- 
rithm that efficiently obtains values of the Hurwitz zeta in the critical strip s = a + it with 
0<O< 1. 

The Hurwitz zeta function may be defined as 



«"»-£ OTPS? 

This series may be directly summed for a > 1, although, for high-precision work, conver- 
gence is annoyingly slow unless one has at least a > 2. A globally convergent series was 
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given by Helmut Hasse in 1930[ 10]: 

Although this series is technically convergent everywhere, in practice, the convergence 
is abysmally slow on the critical strip. A not unreasonable algorithm may be given by 
considering a Taylor's series in q. Expanding about q = 0, one obtains 

(i.3) «^)=i+£(-9)"f i+ r 1 )^ +w) 

It is not hard to show that the sum on the right is convergent for \q\ < 1. Alternately, an 
expansion may be made about q = 1/2: 

c ('. * + = £(-*)" ( 1 + r 1 ) ( 2s+n - o «* + ») 

It may be shown that the above has a region of convergence \q\ < 1 /2. Either of these ex- 
pansions provide a computable expression for the Hurwitz zeta function that is convergent 
on the entire complex s-plane (minus, of course, s = 1, and taking the appropriate limit, 
via l'Hopital's rule, when s is a non-positive integer). The principal drawback to the use of 
these sums for high-precision calculations is the need for many high-precision evaluations 
of the Riemann zeta; the compute time can grow exponentially as the number of required 
digits is increased, especially when q approaches the radius of convergence of the sums. 

It is the poor performance of these sums that motivates the development of this paper. 
Since the Generalized Riemann Hypothesis can be phrased in terms of the values of the 
Hurwitz zeta function on the critical strip, it is of some interest to have a fast algorithm for 
computing high-precision values of this function in this region. There seems to be a paucity 
of work in this area. There is a discussion of fast algorithms for Dirichlet L-functions in 
|[T6l (the author has been unable to obtain a copy of this manuscript). 

The Hurwitz zeta may be expressed in terms of the polylogarithm| 12| (sometimes called 
the "fractional polylogarithm" to emphasize that s is not an integer): 

(1-4) LU(z)=l4 

n=\ 

by means of Jonquiere's identitv lfl2l Section 7.12.21 lfTTl 

(1.5) Li. (J**) + (-iyLi s (e- 2 ™«) = -s,q) 

and so one might search for good algorithms for polylogarithms. Wood[20] provides a 
extensive review of the means of computing the polylogarithm, but limits himself to real 
values of s. CrandallJT) discusses evaluating the polylogarithm on the entire complex z- 
plane, but limits himself to integer values of s. Thus, there appears to be a similar paucity 
of general algorithms for the polylogarithm as well. The series defining the polylogarithm 
may be directly evaluated when \z\ < 1, although direct evaluation becomes quite slow 
when \z\ > 0.95 and a < 2. There do not seem to be any published algorithms that may be 
applied on the critical strip. 

The primary effort of this paper is to take the algorithm given by Borwein||4][5), which 
is a simple Pade-approximant type algorithm, and generalize it to the polylogarithm. The 
result is a relatively small finite sum that approximates the polylogarithm, and whose error, 
or difference from the exact value, can be precisely characterized and bounded. Increased 
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precision is easily obtainable by evaluating a slightly larger sum; one may obtain roughly 
2N to AN bits of precision by evaluating a sum of N terms. The sum may be readily 
evaluated for just about any value s on the complex s-plane. However, it is convergent only 
in a limited region of the z-plane, and specifically, is only convergent when 



This is sufficient for evaluating the Hurwitz zeta for general complex s and real < q < 

I . Unfortunately, there does not appear to be any simple and efficient way of extending 
this result to the general complex z-plane, at least not for general values of s. By using 
duplication and reflection formulas, one may extend the algorithm to the entire complex 
z-plane; however, this requires the independent evaluation of Hurwitz zeta. 

Although the Hurwitz zeta function may be computed by evaluating the Taylor's series 

I I. 31 directly, a superior approach is to perform an Euler-Maclaurin summation (thanks to 
Oleksandr Pavlyk for pointing this out lfMl ). The summation uses the standard, textbook 
Euler-Maclaurin formula|fl] 25,4,7], and is applied to the function f(x) — (x + q)~ s . How- 
ever, the application is "backwards" from the usual sense: the integral of f(x) is known (it 
is easily evaluated analytically), and it is the series, which gives the Hurwitz zeta, which 
is unknown. All of the other parts of Euler-Maclaurin formula are easily evaluated. The 
result is an algorithm that is particularly rapid for the Hurwitz zeta function. It outperforms 
direct evaluation of the Taylor's series by orders of magnitude. It is faster then evaluating 
ll.5l (which can be computed for real values of q). 

The development of this paper proceeds as follows. First, a certain specific integral rep- 
resentation is given for the polylogarithm. This representation is such that a certain trick, 
here called "the Borwein trick", may be applied, to express the integral as a polynomial 
plus a small error term. This is followed by a sketch of a generalization of the trick to the 
convergence acceleration of general oscillatory series. The next section selects a specific 
polynomial, and the error term of the resulting approximation is precisely bounded. This 
is followed by a very short review of the application of the Euler Maclaurin summation. 
This is followed by a brief review of the duplication formula for the Hurwitz zeta, and a 
short discussion of ways in which this numerical algorithm may be tested for correctness. 
Measurements of the performance of actual implementations of the various numerical al- 
gorithms is provided. This is followed by a detailed derivation of the monodromy group, 
and a discussion of Apostol's "periodic zeta function". 

The algorithm has been implemented by the author using the Gnu Multiple Precision 
arithmetic library [8'|, and is available on request, under the terms of the LGPL license. The 
paper concludes with some intriguing images of the polylogarithm and the Hurwitz zeta 
function on the critical strip. 



2. The Polylogarithm 

The polylogarithm has a non-trivial analytic structure. For fixed s, the principal sheet 
has a branch cut on the positive real axis, for 1 < z < +°°. The location of the branch- 
point at z = 1, as is always the case with branchpoints, is the cause of limited convergence 
of analytic series in the complex z-plane. This is noted here, as this branch point has a 
direct impact on the algorithm, preventing convergence in its vicinity. Besides this, the 
polylogarithm has many other interesting properties, which are not reviewed here. 

The development of the algorithm requires the following integral representation. 
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Lemma 2.1. The polylogarithm has the integral representation 



Ll s (z) = ^r-r / — : dy 



r(i) Jo 1—yz 

Proof. This identity is easily obtained by inserting the integral representation of the Gamma 
function: 

Li *(z) = fe-^du 
F(s) ^ n s Jo 



J-£z" fe^f-'dt 



r(i) 7 
1 



r(s) 7o 



i 



1 



dt 



and then finally substituting t = — logy in the last integral. Although this derivation ca- 
sually interchanges the order of integration and summation, one may appeal to general 
arguments about analytic continuation to argue that the final result is generally valid, pro- 
vided one is careful to navigate about the branch point at z = 1 . □ 

A more sophisticated presentation of this theorem is given by Costin[6] Thm. 1]. 

3. The Borwein trick 

The Borwein trick uses the integral form to find a simple, finite sum that approximates 
the desired value arbitrarily well. The trick consists of two steps. The first is to write the 
above integral in the form 

U s ,z) = 777TTTTT / \ dy 

f{\/z) T{s) Jo 1—yz 

(3-D = ".W^i/'MaW-* 

f(i/z) r(s) Jo 1—yz 

The second step is to find a sequence of polynomials p n (z) to be used for f{z), such that 
the integral on the right hand side can be evaluated as a simple, finite sum, while the left 
hand side can be shown to bounded arbitrarily close to zero as n — > °o. The right hand side 
may be easily evaluated, by employing the following theorem. 

Lemma 3.1. Given a polynomial p n {y) of degree n, it can be shown that 

(3.2) r n (y) = ^ 

1—yz 

is again a polynomial in y, of degree n — 1, for any constant z, provided that z ^ 0. 

Proof. This may be easily proved, term by term, by noting that (x n — a n )/(x — a), with 
a — 1/z, has the desired properties. □ 

An explicit expression for r„ (y) is needed. Write the polynomial as 

n 

Pniy) = 

k=0 
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while for r n assume only a general series: 

(3.3) r n (y) = Y d c k y k 

k=0 

Setting y — 0, one immediately obtains cq = «o — Pn ( 1 /z) ■ Higher coefficients are obtained 
by equating derivatives: 

1 — yz L 

where r„ (y) is the £'th derivative of r„ with respect to y. Setting y = in the above, one 
obtains the recurrence relation = cik + zck-i which is trivially solvable as 



Ck 



= z k 



1\ * a; 



From this, it is easily seen that c„ = and more precisely c„+ m = for all m > 0. This 
confirms the claim that r„ (y) is a polynomial of degree n — 1 in y. 
This is immediately employed in equation l3.1l to obtain 

To obtain a good approximation for Li s {£), one needs to find a polynomial sequence such 
that t, goes to zero as n — > °o for the domain (s,z) of interest. That is, one seeks to make 

1 z f l p„(y)\logy\ s - ] 



(3-5) S( JjZ ) = _^ - / 1 Jy 

Pn(l/z) r(i) Jo 1 - yz 

as small as possible. One possible sequence, based on the Bernoulli process or Gaussian 
distribution, is explored in a subsequent section. 

4. Convergence Acceleration of Oscillatory Series 

For the case of z — —1 and p n (y) — y'\ the technique developed above corresponds to 
Euler's method for the convergence acceleration of an alternating series[l, eqn 3.6.27]; 
a particularly efficient algorithm for Euler's method is given by van Wiingaarden |fl9l . 
Polynomials that provide much more powerful series acceleration are discussed by Cohen 
efa/|[T8 1, but again in the context of z = — 1 ; these include the Tchebysheff orthogonal poly- 
nomials considered by Borwein, which do a rather good job of minimizing the integrand of 
equation l3.5l when z = — 1. In a certain sense, an alternating series can be considered to be 
a series from which an explicit factor of z" has been factored out, for the very special case 
of z = — 1 . Thus, one is lead to ask if one can find series acceleration methods for more 
general oscillatory series, where the oscillatory component can be thought of as going as 
z n for some root of unity z- 

The Euler method applied to the Hurwitz zeta function essentially leads to the Hasse 
series 1 1.2l This may be seen by noting that the inner sum is just a forward difference: 



!(-!)*( I )(q + ky- s = A»q l - 



where A is the forward difference operator. The Hasse series has improved convergence 
in that the result is globally convergent, whereas the traditional series representation of 
equation 11.11 is not. However, a quick numerical experiment will show that the Hasse 
series does not converge rapidly, especially in the critical strip. 
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The poor convergence of these traditional techniques motivates the development of this 
paper. Strictly speaking, the slow convergence can be attributed to the fact that the naive 
series representations [TTT1 or 1 1 .41 are not strictly alternating series, as would be required to 
justify the application of the Euler method. Instead, the series are a superposition of os- 
cillatory terms; for the polylogarithm, an explicit oscillation coming from z" = \z\ n e mmgz 
and a slower oscillation coming from vT s = «~ a e~ !Tlog ' 1 . Thus, what one needs is a theory 
of sequence acceleration not for alternating series, but for oscillatory series. The poly- 
logarithm approximation represents the ad-hoc development of a special case for such a 
general theory. 

A loose sketch of a possible general theory is presented below. Suppose that c„ is 
some more or less arbitrary sequence of numbers, oscillatory in n, and that one wishes to 
compute the sum 

S= E c " 

n=0 

but that one wants to apply convergence acceleration techniques to its evaluation. If the 
leading component of oscillation is c„ ~ cosnG, then one makes the Ansatz that there must 
exist some dual series s„ oscillating as s„ ~ sinnO to give e„ = c„ + is„ ~ e m% . One may 
then consider a sum of a„ = e m ^^e„ and consider the summation of an alternating series 

I(-iy«« 

n=0 

By construction, the a„ are now presumed to be relatively tame, so that series acceleration 
techniques can be applied. This trick may be applied to any oscillatory sequence; one need 
not know a precise or a priori value for 0; one need only guess at a reasonable value, based 
on the data at hand. 

Alternately, one may consider replacing the forward difference operator Ab n — b n+ \ — 
b„ by A q b„ = b n+ i — qb„. Straightforward manipulations lead to 

oc / yn+1 
m=0 V 1 1 Z / n=0 

with the traditional Euler's series regained by taking q = 1 and z = — 1. Assuming that 
the oscillatory part z is known, or can be approximately guessed at from the period of the 
oscillation, then choosing q = — l/z should make the left hand side an accelerated series 
for the right hand side. However, there are more powerful techniques. 
Suppose one is interested in the sum 

(4.1) S(z) = £z"fc„ 

n=Q 

In general, one is interested in general sequences b„. Suppose, however, that the b„ are 
expressible as moments, so that 

bn= f 1 y"g(y)dy 
Jo 

for some function g(y). Since g(y) is still fairly general (subject to constraints discussed 
below), this assumption does not overly restrict the b„. If the b„ are expressible as mo- 
ments, then one has 

Jo l—yz 
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Applying the Borwein trick, one obtains 



with 



and 



S„(z) = 



S(z)=S n (z)H{z) 

" l Pn{y)~Pn{\) 



l r l p,Ay) 

Pn (i) Jo 1 



-yz 



Uz) = 



Pn(y)g{y) 
Pn( l z )Jo l-yz 



g(y)dy 



dy 



The goal is to show that the S„ (z) are easily evaluated, while also showing that the £,„ (z) 
are arbitrarily small. The first is easy enough: using equations [32] and [33] one has 



S„(z) 



1 



n-1 



Y, C k b k 



With a suitable choice of p n {z), the coefficients Ck are presumably easy to evaluate. To 
show that the S„ (z) is really a series acceleration for the geometric series 14.11 one must 
show that |^„(z)| is bounded. If g(y) > 0, then 



\Uz)\ < 



Pn (yo) 



Pn < 



\S(z)\ 



where yo is the point on the unit interval where p n (y) assumes its maximum. Provided that 
one can find a polynomial sequence such that 



Pn (yo) 



Mi) 

for some number A > \l/z\, then one has that the series S„(z) converges to S(z) more 
rapidly than the simple partial sums K = nZ*fr* of equation |4. 1 1 

Note that the above is a proof of a series acceleration method for a more-or-less general 
sequence of b n , and is in no way particular to the polylogarithm or the Hurwitz zeta. The 
only ingredients of the proof were that the b n are "well-behaved" - in this case, being 
expressible as moments of a rather general function giy). For the proof to hold, g(y) 
needs to be (mostly) positive, and integrable; it need not be analytic, differentiable or 
even continuous. It needs to be "mostly" positive only so that the bound on the integrand, 
Pn{y)g{y) < \Pn (yo)\g(y) can be established. Thus, one might hope that the acceleration 
method might work, even for those general cases where g(y) is not explicitly known, but 
the b n are somehow "reasonable". 

For the special case of z — — 1 , i.e. for the case of an alternating series, Cohen et al. [ 1 8 1 
suggest some remarkably strongly converging polynomials, with values of A from 5.8 to as 
much as 17.9. The generalization of those results to arbitrary z is not immediately apparent, 
but is surely possible. 

There is another, possibly fruitful direction one might take. This hinges on noting that 
the integrands above take the general form of a Pade approximation, and one proceeds 
from there. 



5. The Gaussian Distribution 



Returning to the polylogarithm, the task at hand is to find a polynomial sequence that 
is suitable for establishing a tight bound on the error term A good bound minimizes the 
integrand, while also maximizing the value of \p,,(l /z)|. In order to suppress the logarith- 
mic branch point at y = in the integrand of 13.51 this section will consider the polynomial 
sequence p„(y) = y"(l —y)"~ a , a sequence which goes over to the Gaussian distribution 
for large n. 

Consider first the polynomial sequence P2n(y) —y n {^ — y)"- It becomes sharply peaked 
at y = 1/2, and subtends a diminishingly small area. For large values of n, this polynomial 
sequence approximates a Gaussian: 
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P2n(y) = ^- 4 ^'- 1/2)2 




1 



Using this in the integral 1331 and assuming that the real part of z is not positive, it is not 
hard to deduce the very crude estimate 



z 



,4(2-1), 

which confirms that t,(s,z) does indeed get suitably small for a certain region in the com- 



plex z-plane. However, in order for equation 13.41 to be useful computationally, an upper 
bound on the value of t, needs to be given, as a function of z and n. This bound is derived 
in the next section. 

The polynomial coefficients are easily obtained, and are 



Cn+k 



= (-If 



a„-i = 

for < k < n 



This leads to 



,2«+l 2n-l 



Li s (z) 

where the c\ are given by 



(z-iy to ( k + l ) s 



z-l 



k—n 



.7=0 



where the summation above is to be understood to be zero when k < n. The above summa- 
tions can be re-organized into the more suggestive form 



(5.1) 



Li s (z) 



7 k 

E- 



1 



2/7 



-k 2n—k 



k=l ' 



a-*r kt n 



The error term t,{s,z) is negligible for only a very specific area of the complex z-plane. The 
region where t, may be ignored as a small error is derived in the next section. Substituting 
z = — 1 in the above formulas agrees with expressions previously given by Borwein[4]. 

The above formula has been implemented using the Gnu Multiple Precision arithmetic 
library (GMP)|[8l, and has been numerically validated for correctness in several different 
ways. The source code, under the license terms of the LGPL, may be obtained by contact- 
ing the author. 



6. Bound on the Error Term 



In order for equation |5. II to be useful computationally, an upper bound on the value of 
\ as a function of n needs to be given. To compute the polylogarithm to some desired 
precision, one infers a suitable value of n based on this bound. However, to achieve this 
desired precision, one must not only choose n large enough, but one must also maintain 
a somewhat larger number of digits in the intermediate terms, as the appearance of the 
binomial coefficient in the equation |5 . 1 1 implies that intermediate terms can become quite 
large, even while mostly cancelling. This section derives an upper bound on the size of \, 
and briefly discusses the issue of the required precision in intermediate terms. 

The general behavior of the integrand appearing in equation [3 . 5 I depends on the whether 
3is > 1 or not; estimates are presented for these two cases. Writing s = a + ix for real a 
and x, and assuming o > 1 and choosing n so that a <n, one has 



P2n(y)\logy 



1-1 



< P2n(y) \logy 



a- 1 



\y(l-y)logyr l (y(l-y)) 



\n-o+l 



Each part of the right hand side is bounded by a Gaussian, and since the product of Gaus- 
sians is a Gaussian, so is the entire expression. From the Stirling approximation, one has 
the well-known identity 



(y(l-y)) a <iexp-4a(y-~ 



The other part is also bounded by a Gaussian 



b(l -y)logy| a < bo(l -yo)logy | a exp 



a(y-yo) 
2yo(l -2yo) 



which is centered on the maximum of |yo (1 — yo)logyo|, that is, at yo = 0.23561058253 . . ., 
where yo is the solution to 

= -f y(l-y) |logy| 
dy 

Curiously, the numerical value of the factor in the Gaussian is 1 /2yo(l — 2yo) = 4.013295587 . 

To evaluate the integral in 13.51 one also needs a bound on the denominator. This is 
furnished by 



l—yz 



<C(z) = 




if5Rz<0 

if0<9k and \z\ < 1 
if0<9tz and Id > 1 



Thus, the integrand is bounded by a Gaussian. Multiplying the two Gaussians, completing 
the square, and evaluating the resulting integral is a bit tedious. The result is that 



L 



P2«(y)|logyr 
l—yz 



- dy 



<C(z) 
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|4yp(l -yo)lo gy | 
4« 



o-l 



■G(s) 



The constant may be taken as \4-yo (1 — yo)logyo| = 1.041381965 . . .. The factor G(s) is 
bounded by 

8yo(l-2yo) 



G(s) < J£exp-(l-2y ) 2 (o-l /c ^ ^ G , 



8yo(l-2y ) + ^ 



+ 1 



< 
< 1 



— exp-(l-2y ) — — 
An a 



when 1 < o. Useful for estimation is that 8yo(l —2y ) = 0.996687115. 
these results, one obtains 



Combining 



(6.1) 



l^,z)l< 



z-1 



C(z)G(s) (1.0414...) 



o-l 



for the case where 1 < a and n chosen so that a <n. 

The case where a < 1 may be evaluated as follows. One writes 



P2n(y)\logy\*~ 



la— l 



< P2n(y) |logy|' 

' V(l-30" +B - 1 



logy 

< /(l-y)"* 0-1 

The integrand may now be recognized as the Beta function, so that 



P2n{y)\logy 
l — yz 



s-1 



■ dy 



<C(z) 



r(» + i)r(»+q) 

T(2« + a+l) 



In the region where \a\ <C n, one may approximate 



l-a 



r(» + i)r(»+a) 2 

r(2n + a+l) ~ 4" 

Combining the various pieces, one obtains a result remarkably similar to the previous one; 
namely, when a < 1 but n is such that \a\ <C «, one has 



(6.2) 



z-l 



)l-a. 



t(z) 



i-l 



For evaluations on the critical line a — 1/2, one needs a good estimate for |r(j) 
which can become quite large as the imaginary part of s increases. A useful estimate for 
|r(j)|~ is given by Borwein EH), for the case where a > —m +1/2: 



1 



1 



|rWI 



< 



i 



n > 



o 



(a+kY 



sinhTTi 

TJt 



n i 



k=0 



(*+*)' 



A remarkable side-effect of the estimation is that, for s equal to negative integers, one 
has that F(s) is infinite, thus implying that the error term is zero. This somewhat surprising 
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result corresponds to the fact that Li_„ (z) has an exact expression as a ratio of two poly- 
nomials, the denominator being of degree n + 1 . Indeed, setting 5 = and so n = 1 into 
equation |5. II gives the exact result 

Li (z) = -i- 
l — z 

while for s = — 1, one must use n = 2, to obtain the exact result 

(1 -z) 

The generic form for the polylogarithm at the negative integers is 

where / ^ X denotes the Stirling numbers of the second kind. 

In general, the formula 15.11 seems to have the curious property of always evaluating 
to exactly the same rational function whenever s is a non-positive integer, and — s < n. 
That is, for fixed non-positive integer s, the sum is independent of n, provided that n is 
large enough. Thus, this polynomial approximation has the pleasing result of giving exact 
answers in the case when an exact polynomial answer is possible. 

For the general case, one may conclude that the estimate [5. ll can be effectively applied 
whenever 
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To obtain a value of the polylogarithm to within some given numerical precision, one must 
invert the formulas 16 . 1 1 or |6T2l solving for the value of n which is to be used in equation 
15.11 To obtain a fixed number of digits of precision, one must carry out intermediate 
calculations with at least that many digits of precision. In fact, one must have more. The 

appearance of the binomial coefficient in equation |5.1| is the problem. Since 2" > 

one concludes that n additional binary bits of precision are needed with which to carry out 
the calculations. As a result, it becomes difficult to implement the algorithm robustly using 
only double-precision arithmetic . For values of z near z = — 1, the algorithm is well-enough 
behaved; however, round-off errors significantly disturb the calculations as z approaches 
+1. This trouble with precision can be at least partly alleviated by making use of the 
duplication formula, as discussed in the next section. 




7. Euler-Maclaurin Summation 



An alternative algorithm for computing the Hurwitz zeta may be based upon the Euler- 
Maclaurin series. This algorithm proves to be quite rapid and efficient|l4). The Euler- 
Maclaurin series may be written as 



if(k) = £/(*) 

k=0 k=0 
P 



(7.1) 



AN) 



f(x)dx 



2k+l 



- (2k)\ dx 2k + 



-fix) 



R 



x=N 



This formula is simply applied to the function f(x) = (x + q) s to obtain the Hurwitz zeta. 
The derivative is particularly easy to evaluate: 

d 2k+l 1 _ _s(s+l)(s + 2)---(s + 2k) 
dx 2k+l (x + q) s ~ (x + q) s+2k+l 

and the integral is trivial: 

-dx = 



In (x + q) s ~" (s-l) (N + qy- 1 
The error made by this expansion is embodied in the term R. It is directly given by 

R=T i_r f i P ^ [x) B ^ x -^ dx 



where the Bk(x) are the Bernoulli polynomials. There are various formulas in the literature 
which may be readily applied to bound the size of this term(T)- Thus, for example, one has 

R< 2 \s(s+ !)■■■(* + 2p)\ 1 
" (2k) 2 p %s + 2p (N + q)* s + 2 P 

Thus, the evaluation of the sum l7.1l merelv requires a suitable choice of N and p to be used. 
As it happens, the situation is even simpler than this. The sum |7.1| is an asymptotic series, 
and it is well-known that the best approximation for an asymptotic series occurs when one 
stops the summation at the smallest term in the series. Thus, it is sufficient to choose a 
value of N, while p is found dynamically as the algorithm runs. 

There remains the question of what value of N to use. This remains an open problem, 
which will not be explored here. However, it seems that an adequate heuristic for most 
common cases is to choose N = D/2+ 10 where D is the desired number of decimal 
digits of precision. The performance of this algorithm is compared to that of the Borwein 
algorithm in the next section. 

XXX to do: characterize the region of convergence for this algo. 

8. Performance 

The performance of both the Euler-Maclaurin summation and the Borwein algorithm 
appears to be very good. As a general rule, both have better performance than the di- 
rect summation formula [T~4l although this depends on whether the poly logarithm or the 
Hurwitz zeta is the desired quantity. It also depends on whether repeated evaluations are 
being done for fixed s but varying q or z, as holding s fixed allows the results of intermedi- 
ate calculations to be cached. The figures [Ol and W2\ compare the performance of actual 
implementations . 

The first figure shows "cold cache" performance, measured in seconds, as compared to 
the number of desired decimal places of precision. The figure is termed "cold cache", in 
recognition of the fact that some constants may be pre-computed. For example, if s is held 
fixed, while z is varied, then the values of n~ s appearing in both the direct sum and the 
Borwein algorithm may be computed once, and then re-used for subsequent calculations. 
As the figures indicate, computing rC s can be very expensive for general complex-valued 
s, and so the caching strategy offers a big gain when $ is held constant. 

As the "cold cache" figure demonstrates, the Borwein algorithm is faster than direct 
summation. The problem with direct summation is that it requires more values of rC s to 
be computed to achieve comparable precision. 

As the "warm cache" figure shows, the Borwein algorithm will still win against direct 
summation, even when the values of the rC s are pre-computed. Even when these are 
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Figure 8.1. Polylogarithm cold cache performance 
Compute time for Li_{0.5+i1 4.1 3} (0.4+i0.3) for cold cache 
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This figure shows the time, in seconds, needed to compute the polylogarithm to the 
indicated number of decimal places of precision. The values chosen for evaluation are 
s = 0.5 + / 14. 134725, which is very near to the first of the non-trivial Riemann zeroes, 

and z = 0.4 + z"0.3, which is relatively near to the origin. For such a small value of z 
(\z\ = 0.5), one might hope that direct summation might compete favorably against the 
other two algorithms; and in particular against the Borwein algorithm, as the point is 
closer to the troublesome branch point at z — + 1 than algorithmically optimum z = — 1 . 
Nevertheless, |z 2 / (z — 1) h=s 0.373 for this point, and the Borwein algorithm wins. The 
Euler-Maclaurin algorithm is relatively disadvantaged, as it needs to be evaluated twice, 

and a value for T(s) must be computed as well. 
The primary contribution to the calculations is the evaluation of rC s in the summations. If 
these are pre-computed, the Borwein algorithm still wins, as shown in the next figure. All 

calculations were performed on a contemporary desktop computer, using the GMP 
library; the algorithms were implemented in the C programming language. Note that the 
actual performance depends on the underlying implementation of log, exp, sin, atan, sqrt, 
gamma and the like; the actual implementation used here is fairly naive and untuned. 



pre-computed, and can be pulled from the cache, direct summation still requires more 
operations. 

A different but equally dramatic set of results obtain, when one compares the perfor- 
mance of the Borwein algorithm applied to the Hurwitz zeta function, as compared to the 
use of the Taylor's expansion ll.3l for the Hurwitz zeta. In this case, as the figures show, the 
Taylor's series outperforms the Borwein algorithm by a constant factor of three, when both 
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Figure 8.2. Polylogarithm warm cache 
Compute time for Li {0.5+i1 4.1 3} (0.4+i0.3) for warm cache 
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This graph compares the compute times for the Borwein algorithm, the direct summation 
of the polylogarithm, and the Euler-Maclaurin series, for a "warm cache". That is, an 
array of values of rC s were pre-computed, prior to beginning the timing measurements. 
These pre-computed values were used in the Borwein and direct sums. For the 
Euler-Maclaurin sum, the value of F(s)/(2kY is pre-computed and cached as well. 
Despite the use of such cached, pre-computed values, the Borwein algorithm still wins 
over direct summation, and appears to be tied with the Euler-Maclaurin sum. This graph 
uses the same s and z values as in the previous graph. 



algorithms use pre-computed constants. However, the cost of pre-computing these con- 
stants skyrockets for the Taylor's series, making it unattractive as the number of required 
decimal places increases. 

Although the Borwein algorithm seems to show a slight advantage over the Euler- 
Maclaurin series when it is used to evaluate the polylogarithm, much of that advantage 
disappears when evaluating near the branch point z = 1 . The Borwein algorithm breaks 
down near the branch point, and requires the (possibly recursive) application of the dupli- 
cation formula (given in the next section) to obtains points sufficiently far away from z — 1 . 
By contrast, the Euler-Maclaurin series seems to happily converge in this area, and so no 
extra steps are required. 

In conclusion, it appears that the Euler-Maclaurin formula provides the best algorithm 
for evaluating the Hurwitz zeta, and can often tie and sometimes outperform the Borwein 
algorithm for the polylogarithm. 
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Figure 8.3. Hurwitz zeta cold cache performance 
Compute time for zeta (0.5+H4.13, 0.2) for a cold cache 
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This figure compares the performance of the Borwein algorithm, the evaluation of the 

Taylor expansion, and Euler-Maclaurin summation, for the Hurwitz zeta. The vales 
computed are for s = 0.5 + i 14. 13 and q — 0.2. For such a small value of q, one might 
have hoped that the Taylor expansion might converge quickly. This appears to not be the 
case, as the binomial coefficients can grow to be quite large, and thus force many terms to 

be summed. 

The figure is for a "cold cache", assuming that no constants have been pre-computed. The 
lines are not parallel, with the Taylor expansion getting progressively worse. Both the 
Borwein and the Euler-Maclaurin algorithms are orders of magnitude faster than the 
Taylor's series for high precisions; the Euler-Maclaurin algorithm appears to be an easy 

winner for all precisions. 



9. Duplication formula 

The region of applicability of the algorithm may be extended by making use of the 
duplication formula for the polylogarithm or periodic zeta. The duplication formula, or 
more generally, the multiplication theorem, is an extension of the well-known Legendre 
duplication formula for the Gamma function[U (6.1.18),(6.3.8)(23.1.10)], into the domain 
of the polylogarithm. 

Thus, for example, the formula lBTTl to gether with the error bounds l6. U6.2l allow F(q;s) = 
Li s (e 2mq ) to be computed for real q in region l/4<^<3/4. To obtain values on the re- 
gion < q < 1/4, one applies the duplication formula 

F(q; S )=2 1 - s F(2q;s)-F(q + ±;^ 
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Figure 8.4. Hurwitz zeta warm cache performance 
Compute time for zeta (0.5+i1 4.1 3, 0.2) for a warm cache 
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This figure compares the evaluation times for Euler-Maclaurin summation, the Borwein 
algorithm, and the Hurwitz zeta Taylor's expansion, for the "warm cache" scenario. In 
this case, it is recognized that if s is assumed to be held fixed, then the values of n~ s 
appearing in the Borwein algorithm may be pre-computed. Similarly, the values of 

/ s-\-n — 1 \ 

£(s + n) and the binomial coefficients I ' I appearing in the Taylor's series and 

the Euler-Maclaurin formula may be pre-computed. 
For this value of q, the Borwein algorithm requires four evaluations of the polylogarithm; 
values of q closer to 0.5 would require only two. The Taylor's series evaluation appears to 
be about eight times faster than these four evaluations (or four times faster than the 
minimum of two evaluations). 
Note that the time axis here is in milliseconds, not seconds. Comparing to the previous 
graph, it is clear that performing the pre-computations can be terribly expensive. The first 
evaluation of the function can take 100 or 1000 times longer than subsequent evaluations 

at the same value of s. 



recursively until one obtains values of q > 1/4. Rearranging terms, one obtains a similar 
formula for iterating values of q > 3/4 until they are less than or equal to 3/4. Thus, 
as q approaches or 1, the algorithm requires more time, but only logarithmically so, as 
— log 2 g or — log 2 (l — q) as the case may be. 

The duplication formula follows from a general p-adic identity 

m=0 V " / 
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which is valid for any positive integer p. One of many ways of obtaining this formula is by 
noting that the function F (s;q) is an eigenvalue of the /?-adic Bernoulli transfer operator 
associated with eigenvalue y 11171 . 

The equivalent formulas! 12, Sections 7.3.1, 7.12.1] for the polylogarithm are 

Li s (z)+Li J (- Z )-2 1 -' s Li J (z 2 ) 
whereas the p-adic version has the appearance of a Gauss sum: 

£ Li s (zt***/') =p 1 - s Li s ( z P) 

m=0 V ' 

The application of the duplication formula, together with the inversion relation [T31 can 
be used to extend the evaluation of the polylogarithm to the entire complex plane. For 
numerical work, both formulas must be applied, one alone is not enough. Consider first 
the application of the duplication formula only. It may be used to take points that are 
near to z = 1, and map them further away from z = 1, into the convergent region for the 
Borwein polynomial. Repeated application allows arbitrarily close approach to z = 1 from 
the left-hand side. The resulting region of convergence is kidney-shaped, with the cusp of 
the kidney at z = 1 , and the kidney containing the unit disk \z\ < 1 . The precise shape of the 
kidney depends on the number of terms one wishes to use in the polynomial approximation. 
The shape that can be achieved while still maintaining good running time is shown in figure 
l~4l However, as can be seen, this strategy barely penetrates the 3iz > 1 region. 

To extend the algorithm to the entire complex z-plane, one must make use of the Jon- 
quiere's inversion formula! 15, pp 27-32] 

(9.1) e-^U s (z)+e^U s (-)=^(ls, '° gZ 



r(s) 3 V 2ni 

together with some sort of independent means of evaluating the Hurwitz zeta function. The 
Taylor expansion 1 1.3 1 is particularly well-suited, as it is rapidly convergent in the vicinity 
of z = 1. Specifically, it is convergent when |logz| < 2tc, although the region of acceptable 
numerical convergence is smaller, roughly |logz| < 71. Either way, this encompasses a 
rather large region in the vicinity of z = 1, which is exactly the region where the Borwein 
algorithm has trouble. The inversion formula is used to move points \z\ > 1 from outside 
of the unit circle, to the inside. Not all points inside the unit circle are directly accessible 
to the Borwein algorithm; the duplication formula must still be used to handle points in 
the disk interior that are near z = 1. Similarly, points with very large z, such as those for 
which |z| > e 2lt , require the use of the duplication formula to be moved into the region 
of convergence for the Hurwitz Taylor series. A typical view of the polylogarithm on the 
complex z-plane is shown in figure [31 

As a practical computational matter, rather than using 19.11 it seems to be more conve- 
nient to use the formula 



T(s) 

(9.2) Li!_ s (z)- 



Although this is equivalent to the equation 19.11 immediately above, it correctly captures 
the appropriate branch of the logarithm that should be used: equation |9. 21 works well for 
values of s in both the upper and lower half-planes, whereas equation |9.1| is more difficult 
to correctly evaluate in the lower half s-plane. In either case, one must use the logarithm 
with an unusual branch cut, taking it to run from z — to the right along the positive real 
axis, as opposed to the usual cut taken for the logarithm. 
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It is of some curiosity to note that the sheet that is forced takes the form 

which is also noted by Cosfin||6] eqn 14] as a Mittag-Leffler type decomposition. A general 
discussion of the monodromy is given in a later section. 



10. Testing and Validation 

The correctness of a given numerical implementation can be validated in a number of 
ways. For non-positive integer s, one has the exact rational functions previously mentioned. 
For positive integer n, one has the relationship 

Li„ (e 2 ^) + (-l)"Li„ (e- 2 ^) = ~^B n (q) 

where the B n (x) are the Bernoulli polynomials. For z = — 1, one regains the Riemann zeta 
function 

M.(-i) = 2i=Jn«*) 

For |z| < 1, the defining series [T~4"l is explicitly convergent, and may be directly summed, 
thus offering a fourth check of the correctness of an implementation. Finally, the multi- 
plication theorem can be used to check for consistency. Each of these quantities can be 
computed by independent algorithms, and thus be used to validate the correctness of a 
polylogarithm implementation. 



11. Branch Points and Monodromy 

The principal sheet of the polylogarithm has a branch point at z = 1, and by convention, 
a branch cut is placed along the positive real z-axis, extending to the right. As one moves 
off of the principal sheet, one discovers that there is another branch point at z = 0. The 
resulting monodromy group is generated by two elements, acting on the covering space of 
the bouquet S l V S 1 of homotopy classes of loops in C \ {0, 1} passing around the branch 
points z = or z = 1 . The author is not aware of any simple published discussion of the 
monodromy for the polylogarithm; a much more abstract discussion is given in [6, 9 IT3l[3ll. 
Thus, this section provides a discussion; note that the correct manipulations to move from 
one sheet to another can be somewhat treacherous and confusing. The 

The inversion formula 19.11 suggests a way to move around the branch point at z = 1 . 
Suppose one starts at the point z = x + ie with x real, positive, and greater than one, and 
e some arbitrarily small positive real number; thus z is very near the branch cut of the 
principal sheet. One wishes to compare this to the straddling value at x — ie. Applying the 
inversion formula, one can bounce these two points inside the unit circle, where there is no 
cut, and thus the polylog differs by (e). The difference across the cut is thus 

A = Li. s (x + ie) — Li j (x — ie) 

;../, (2tc)' s L A lnx + /e\ „/ lnx-/e\l . . 

Now, since x > 1, a naive application of this formula yields A = since log(x + ie) — 
\og(x — ie) = (e); but clearly A cannot be zero. To resolve this situation, one must make 
the Ansatz that the cut of the logarithm should be crossed. This may be done by taking the 
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Figure 11.1. Polylogarithm Monodromy 
-» * 




A graph of the polylogarithm monodromy, given by equation ll 1.41 



cut of the logarithm to extend to the right, instead of to the left, as it is by convention. One 
then gets 

(11.1) log(x + /£)-log(x-/e) =2niN+o(e) 

for some integer N. The difference across the cut, for positive N, is then 

1 



C(*,«)-CMr- 



AT— 1 



with q = ]n(x)/2%i. By swinging the cut of the logarithm so that it extends to the right, 
one obtains that the real part of q is positive; the real part of q runs between and 1 . As 
a result, the value of (k + q) s is unambiguous; as otherwise taking something to the power 
s also begs the question of which sheet the power must be taken on. That is, for general 
complex w, the power function is multi-sheeted: 

2) yy s pSlnw ^ ^.v(lnvv+27c/M) „2nisM^ 



W 



for some integer M representing the sheet of the power function. Since the real part of q is 
taken as positive, one can temporarily operate with the assumption that M = 0. 

Taking N = 1, the above reasoning provides an excellent description, which may be 
verified numerically. The concentric semi-circles visible in the image [3] can be entirely 



explained by the behavior of q s as q — i 
the N = 1 sheet and the N = sheet is 



0, that is, as z — » 1 . Thus, the difference between 



(113) 



ins/2 



(2< 

ro) 



lnz 
2~Ki 



Properly speaking, Ai is a function of s and z, and so one should write A\(s;z) to signify 
this. However, for ease of notation, this marking is dropped, and is taken implicitly in 
what follows. The difference Li s (z) — Ai is illustrated in figure [7J the concentric rings 
are seen to be cancelled out precisely, leaving behind a smoothly varying function, having 
the expected smooth structure. Moving across the cut, for 3iz > 1, the joint appears to be 



smooth. For general positive N, it is easy to confirm numerically that the discontinuity 
between the A^'th sheet, and the N — 1 'th sheet, across the 3iz > 1 cut is 



The discontinuity across the cut 3iz > 1, for negatived, follows similarly. Next, taking 
N = — 1 in equation 1 1 1.11 (or N = in the equation immediately above, which can serve 
as a point of confusion) the difference is given by 1 + q) = e~ ms / (1 — q) , or 

r(s) V 2jci 

The figure [8] illustrates the difference Li v (z) — A_j. For general negative AT, the differ- 
ence between adjacent sheets is then 

A _ N = e -ins/2(^r iN _ q y-l 

r( s ) 

Again, it can be numerically verified that that the transition from sheet to sheet is smooth 
as one crosses the cut 3iz > 1 . 

It is perhaps more clear to use explicit topological language. Let mi represent the homo- 
topy class of all loops based at some point z on the complex plane, that wind once around 
the branch-point z = 1 in the positive direction. The action of nt\ on the polylogarithm has 
the effect of carrying the polylog from one sheet to the next. In the above discussion, it 
was determined that 

m\ ■ IA S (z) = Li s (z) - Ai 
The logarithm in A] has a branch point at z = 0. That is, after acting with mi , one is now on 
a sheet with a cut extending from z = to the right. Let let mo represent the homotopy class 
of loops that wind once around the branch point z = in a right-handed fashion. Acting on 
the logarithm, one has 

mo Tnz = lnz + 2ni 
Recalling the definition of q = lnz/2ni, one thus has mo ■ q = q + 1, and so 

mo • An — Ajv+i 

The principal sheet of the polylogarithm has no branch point at z = 0, and so one has 

m • Li s (z) = Li,, (z) 
Winding in the opposite direction, one has 



mj 1 • Li,, (z) = Li s (z) - A_ 



l 

In order for mi to be properly considered as the group-theoretic inverse of mi, one must 
have mi • nC[ 1 • Li s (z) = m]~ -mi ■ Li. s (z) = Li, (z). The resolution of this implies mi ■ A_ i = 
— Ai, which in turn implies that m\ ■ q = q+l. This seems strange, as the logarithm has 
no branch point at z = 1, and so there is nothing to wind around. By this argument, one 
would have expected that mi had no effect on q at all. The point is subtle and is worth 
establishing clearly; it is the joining of the polylogarithm to the logarithm cuts that causes 
the effect. Consider only the logarithm (not the polylogarithm), arranged so that the branch 
extends to the right. Consider starting just below the real axis, and winding around in a 
right-handed fashion around the point z — 1, and finishing just above the real axis. There 
is no obstruction at z = 1, and so this loop can be shrunk to a very short line segment. 
None-the-less, the short line segment crosses the cut, and its effect is to move to a different 
sheet. It is for this reason that one has mi ■ q = q+ 1. A different set of group generators 
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will be defined below, that do capture the idea that winding around z = 1 should have no 
effect on q. 

One can now consider the task of more complex paths from sheet to sheet. Passing 
twice around the z = 1 branch, one has 



m\ ■ Li. s (z) = mi ■ [Li s (z) - Ai] = Li,- (z) - Ai - A 2 
Repeating this gluing n times, so that each time, one pastes the sheets so that crossing the 
z > 1 cut is smooth and differentiable, gives 



•Li s (z) 



k=\ 



LL(z) 
LL(z) 

U s {z)-j™ 12 



ins/2 ( 2% T y 

(2nY 



1 



r(s) 



1+9) 
lnz 



l-.v 



271/ 



lnz 

271/ 



which is recognizable from the initial considerations. 
To capture the idea of there being no obstruction at z 



1 for the logarithm, one may 



define a group element g\ =mim , so that one has g\ ■ lnz = lnz. Then, define go = mo. 
In terms of these generators, one has the relations 



go-q 


= q+l 


gi-q 


= q 


go ■ Lis (z) 


= U s (z) 


gi-U s (z) 


= Li s (z) 


go- An 


= A-N+l 


gi -A N 


= 



- Ai 

forJV>0 



To complete the picture, one needs the action of go on A_i, or equivalently g^ 1 on A;. 
There is some ambiguity, and thus, room for confusion, as the term (— l) s arises, which 
may be taken as e ms or as e~ ms , with the last two inequivalent for non-integer s. The 
resolution lies in considering the joining of sheets across the cut that runs between < 
3iz < 1. Visually, the cut can be clearly seen in figure[8] Consider now the task of glueing 
sheets across this cut. For a point z just above the line connecting z = and z = 1, one 
has q = E + iv for some small, positive e and some general, positive v. Just below this line, 
one has q = 1 — e + /v. That is, q differs by 1, of course. This cut cannot be glued to the 
polylog, of course; it must be glued to another sheet of the logarithm. The correct gluing, 
for z—x real, < x < 1, is given by 

J2ns * i(*-fe)] =0 



lim [Ai(jc + /e)+e'' 2m A_ 

e-»0 



That is, 

go-A- l = -e- i2ns A l 

This somewhat strange form is nothing more than an accounting trick; it is the result of 
providing a definition of A_a? that had a "natural" normalization, avoiding a minus sign. 
The price of that definition is this glitch. In all other respects, the homotopy proceeds as 
expected, so that, for example, 

go ■ A_n = A_Af + i 

forJV> 1. 
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The free combinations of powers of the two operators go and g\ (or mo and mi) generate 
a group, the monodromy group of the polylogarithm. For s — ma positive integer, the mon- 
odromy group has a finite-dimensional representation of dimension m+l. Well-known is 
the case for s — 2, the dilogarithm, where the representation forms the discrete Heisenberg 
group. This may be seen as follows. For s = 2, one has 

A„ = 2%i (lnz + (n - l)2%i) 

Repeated applications of go and gi result in a linear combinations of Li2 (z), lnz and 1; 
no terms of any other type appear. One may take these three elements as the basis of a 
three-dimensional vector space. Writing the basis as column vectors, the representation is 



4k e\ = 



— 2;t/lnz i ► (?2 = 



Li 2 (z) e 3 



so that the action may be represented as 





' 1 





" 




' 1 


1 





81 = 





1 


1 


and go = 





1 













1 










1 



These two matrices may be seen to be the generators of the discrete Heisenberg group 
^3(Z). The general group element of the discrete Heisenberg group is 

1 a b 
1 c 
^ 1 

for integers a,b,c. By convention, one writes x for go and y for g\, and defines a new 
element z (having no relation at all to the argument z of the dilogarithm) as the commutator 



z = xyx l y l ,so that 



The element z is in the center, so that it commutes with x and y, so zx - 
These relations may be taken as giving the group presentation, so that 



xz and zy = yz- 



X 3 (Z) = ((x,y)\z = xyx- 1 y 



l-i 



zx = xz. 



zy = yz) 



Curiously, one should note that the discrete Heisenberg group has an alternate presenta- 
tion that is reminiscent of the braid group. This is perhaps not a complete surprise, as a 
monodromy can be generally understood to be a subgroup of the braid group or mapping 
class group of the punctured disk. The mapping class group permutes the punctures of the 
disk in a continuous fashion, so that the result is an intertwining of strands. For a disk 
with two punctures, the braid group is B3 consists of three strands. In this case, two of 
the strands can be visualized as two parallel lines, each representing one of the two branch 
points of the polylogarithm. The third strand is the monodromy path, which weaves about 
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the two straight lines, but must always return to the middle. Now, the braid group Bj has 
two group generators, denoted Oi and 02 by convention, which denote a single twist of the 
left or right pair of strands of the braid. To weave the middle strands of the monodromy 
about the two branch points, one identifies go — > Oj and g\ — > At this point, ignoring 
the polylog, and limiting ones attention to the universal covering space for the homotopy 
group over two branch points, one very simply has that go and g\ are the generators of the 
free group in two letters. Now for the curious resemblance. The braid group B3 has the 
group presentation O1O2O1 = 020102- The discrete Heisenberg group has the remarkably 
similar identity xyyx = yxxy. 
For s = 3, one may write 



1110 
12 
10 
1 



and gi 



10 
10 
11 
1 



for the basis 



4-7t i 1 ^ € \ 



2% lnz 1— ► e2 



-ju'(lnz) 1 ► £3 



Ll3 (z) !-» <?4 = 



The generated group is not the four-dimensional Heisenberg group, as there simply aren't 
enough generators for that. Its not a representation of the three-dimensional Heisenberg 
group, as the commutator w = gogig^g^ 1 is not in the center, since wgo 7^ gow. However, 
one does have wg\ — g\w, and from this one may deduce that the monodromy group has 
the presentation 



(11.4) 



M = (go,gi\gigog 1 g =gog 1 g gi) 



■Mj (w) where 



In any case, the group is solvable and thus thin. 

The Heisenberg group may be regained as the quotient group #3 (Z) 
(w) is the conjugate closure of w: 

(w) = {gwg~ { \geM} 

The conjugate closure is, of course, a normal subgroup of M. Since gi already commuted 
with w, what the quotient group construction does is to impose (w)go = go(w), which 
shows up on the cosets as the other, "missing" relationship needed to get the Heisenberg 
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group. It is thus that the presentation of the Heisenberg group is regained. The general 
form of an element h e (w) may be taken to be 

1 2n + l " 
10 -2 
1 
1 

for integer n. Thus, it is clear that (w) is abelian, and that (w) = Z. 

For general integer s = m, the monodromy group is always unipotent, and thus a nilpo- 
tent group. The generators take the form 

1 ••• " 
1 : 

: '-.0 

1 1 
••• 1 _ 

where C are the binomial coefficients written in an upper-triangular fashion (thus, for 
example, for m — 4, one adds a column of (1,3,3,1) and, for m = 5, another column 
(1,4,6,4, 1), etc.). This form is determined by taking as the basis vectors the monomials 
appearing in the expansion of (lnz + 27t/) m (whence the binomial coefficients), and nor- 
malizing with a factor of 2ni/T(m). The generated monodromy group is isomorphic to the 
above s — 3 case, as it has the same presentation. The element w = gogigQ g± commutes 
with g\, just as before; but gow ^ wgo. Thus, the conjugate closure (w) is again a one- 
dimensional abelian group, that is, (w) = Z, and the resulting quotient group is again the 
Heisenberg group. 

Note that the existence of polylogarithm ladders seems to be equivalent to the statement 
that products of monodromy groups can be factored in interesting and unusual ways. XXX 
Expand on this remark. 

For s not an integer, the action of go does not close, and the general vector-space rep- 
resentation is infinite dimensional. The action of the go monodromy can be understood to 
be the shift operator on an infinite-dimensional vector space, whose basis vectors may be 
taken to be A^; the operator gi acts to inject into this shift space. Defining w as before, it 
is not hard to determine that w ■ An — An and that w ■ 1A S (z) = Li s (z) + Ai — A2. Remark- 
ably, the same commutation relation holds as before, in that g\ w ■ Li. s (z) = wgi ■ Li s (z) — 
Li 4 (z) — A2. This, together with the more trivial result giw-A^ = wg\ ■ An = An shows that 
the infinite-dimensional vector space representation of the monodromy has the same pre- 
sentation as the finite-dimensional case. Thus, one concludes that, for all complex values 
of s ^ 2, the monodromy of the polylogarithm is given by equation l 11.41 

It is not clear whether there may also be smaller, finite-dimensional representations for 
particular values of z. XXX Expand on this. 

The analog of the Dirichlet L-functions are functions of the form 

/^ s Ex(^)Uv(z e 27C,m/ ") 

m= 1 

where %(m) is a primitive character modulo p. These functions have the strange property 
of having branch points at e 2%,m l p whenever %{m) is not zero. The resulting monodromy 
groups have a more complex structure as well. 
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h = go n wg n = 



go 



C 



••• 1 



and gi 



12. The "Periodic" Zeta Function 



The ability to compute the polylogarithm allows one to visualize some of the interest- 
ing sums that occur in number theory. One such example is the so-called "periodic zeta 
function", defined by Apostol[2] Thm 12.6] as 

(12.1) F(q;s) = £ — 

Clearly, one has F(q, s) — Li s (e 2lt,i? ) . The periodic zeta function occurs when one attempts 
to state a reflection formula for the Hurwitz zeta function, as 

(12.2) = [e-^ 2 F(q;s)+e^ 2 F(l-q;s) 

What makes the periodic zeta function interesting is that it is not actually periodic. That 
one might think it is seems obvious from the definition 112.11 the direct substitution of 
q — > q + 1 gives the false suggestion of periodicity in q. This is false because, in fact, 
Li s (z) has a branch point at z = 1 ■ The "periodic" zeta function is multi-sheeted, and, 
attempting to trace a continuous path from q — ► q + 1, while keeping q real carries one 
through the branch-point, and from one sheet to the next. This is illustrated graphically, in 
figure H2J] 

The figure 112.11 shows an oscillatory behavior that clearly diverges as q — > 0. This 
oscillation is can be simply explained. Substituting s = j + ix into equation ll2.2| gives 

C f - - it, qj = ^^^^ \e I F{q;s) +ie / F(l ~q;s) 
after the substitution of 

' 1 \ Jne^ 
- + it) = 

,2 / VcoshTft 

where § is a purely real phase that can be expressed as a somewhat complicated sum| 1 
eqn 6.1.27]. For reasonably large, positive x, such as x = 25 in the picture, one may ignore 
the second term, so that the whole expression simplifies to 

q-l; + iz) =expz(^-<|> + xlog27i) (jQ-fcgJ +0 (C" T ) 

for some real constant C > 0. Then, as q — > 0, the leading contribution to the Hurwitz zeta 
comes from the « = term in equation ll.il that is, q~ s = e n]ogq / y/q, so that 

/ 1 \ e "l°g9 e 'V 

(12.3) F[q;- + n w — as q -> 

V 2 / V? 

for some fixed, real phase \\t that is independent of q. As is eminently clear from both 
the picture [12T1 and from the approximation ! 12.31 the limit of q — * of F(q,s) cannot be 
taken: this is the branch-point of Li s (z) at z = 1 . 

Curiously, the limit q — > 1 does exist, and one has F(l,i) = ^(5) the Riemann zeta. The 
approximation 1 1 2 .31 hints at how to move through the branch point, from one sheet to the 
next: I12.3l is the discontinuity between sheets. That is, one has, for large, positive x 

1 \ / 1 \ e m °^e^ 
q+l; - + ix «F q,-+ii - 



2 J V 2 J ^q 



25 




This graph shows the real and imaginary parts of the periodic zeta function 

F(q, S )=U s (e 2 ™1) 

for s = j + i25. This value of s is close to a zero of the Riemann zeta function, at 
s = j + i25. 01085758014 . . . Thus, both the real and imaginary parts approach zero at 
q = 1, as well as at q = j. The increasing oscillations as q — > are due to the contribution 
of the very first term of the Hurwitz zeta: that is, these oscillations are nothing other than 
that of q~ s = e' 25]ogq / y/q. Subtracting these away is the same as analytically continuing 
to the region q > 1, and matches the coarser set of oscillations, which are given by 
(1 + q)~ s = e a51 °g( 1 +9)/yT+^ 
Noteworthy is that the presumed "periodicity" in q is very misleading: the image suggests 
an essential singularity at q = 0, and continuing, logarithmically slower oscillatory 
behavior for the analytic continuation to the region where q > 1 . 



as the formula that approximates the movement from one sheet to the next. In essence, this 
shows that the "periodic" zeta function is not at all periodic: F{q+\;s)^F(q;s) whenever 
%s< 1. 

The complete relationship between the Hurwitz zeta and the periodic zeta is rather sub- 
tle. For example, if instead one considers large negative x, and graphs the periodic zeta, 
one obtains what is essentially the left-right reversed image of 112. II with oscillations ap- 
proaching a singularity at q — > 1 . One may easily numerically verify that these oscillations 
are precisely of the form (1 — q)~ s . From this, one may deduce that the correct form of the 
reflection formula is 



(12.4) 



,l-s) = 



(2jc)' 



e™^( S ,q)+e-™'X( S ,l-q) 



— ins j '2? 
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which captures the oscillatory behavior at q — > for large positive x and the oscillations at 
q — > 1 for large negative x. 

The constraint of keeping q real causes integer values of q to correspond precisely to 
the branch point of the polylogarithm. This makes reasoning about the continuity of the 
periodic zeta at integral values of q rather turbid. Since the Riemann zeta function lies 
precisely at the branch point, this seems to also make it difficult to gain new insight into 
the Riemann zeta in this way. 

13. Conclusion 

Both the Borwein and the Euler-Maclaurin algorithms appears to offer a stable and fast 
way of computing the Hurwitz zeta function for general complex s and real q, and the 
polylogarithm for general complex s and z. The Euler-Maclaurin algorithm offers superior 
performance for the Hurwitz zeta. An actual implementation using a variable-precision 
library shows that all of the algorithms are quite tractable. 

An unexplored area is the optimal implementation of the algorithms using standard 
IEEE double-precision floating-point mathematics. Preliminary work shows that some of 
the intermediate terms are just large enough so that rounding error may be of concern; 
the region of high-precision convergence has not been characterized. The propagation of 
rounding errors through the algorithm has not been characterized; it may be possible to re- 
arrange terms to minimize the propagation of errors. Such a characterization is necessary 
before reliable double-precision code can be created. By contrast, such considerations can 
be swept under the rug when employing arbitrary-precision arithmetic. 

For the Borwein algorithm, evaluating the bounds |6TI and [6721 so as to obtain an order 
estimate, is non-trivial, and can account for a significant amount of compute time. A 
speedier heuristic is desirable. By contrast, the Euler-Maclaurin algorithm seems to have a 
very simple heuristic that is quite adequate for general-purpose applications. In either case, 
it would be desirable to have a more refined analysis and presentation of the heuristics. 

It may be possible to find even more rapidly converging algorithms, by generalizing 
some of the results from Cohen et a/0~8), or by appealing to the general theory of Pade 
approximants. Both approaches seem promising, although that of Pade approximants may 
yield a more general theory. 

The polylogarithm and Hurwitz zeta functions are both special cases of the Lerch tran- 
scendental 



It should be possible to extend the techniques in this paper to provide a rapid, globally 
convergent algorithm for the Lerch transcendental, thus enabling a deeper numerical ex- 
ploration of its peculiarities. 

This paper concludes with the application of the algorithm to render some pretty images. 
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Appendix: Some Pretty Pictures 

Below follow some interesting pictures of the polylogarithm and the Hurwitz zeta func- 
tion for values of s on the critical line s = 1/2 + i%. They are shown here, as they are not 
commonly known. Each figure represents hours of computation on modern computers. 
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Figure . 1 . The Periodic Zeta 




This image illustrates the so-called "periodic zeta function" 



n=\ ^ 

Graphed is the magnitude of F(q;s), with a color scale such that black represents zero, 

greenish-blue being values on the order of one, yellow on the order of two, with 
increasingly large values shown as orange-red. Along the horizontal axis runs q, taken to 
be real, from q = on the left to q = 1 on the right. Then, writing s = j + ix, the value of 
X is varied along the vertical axis, running from x = at the bottom of the image, to 
X = 50 at the top of the image. Zeroes of the Riemann zeta function =F(l,s) are 
located where the blue lines intersect the right edge of the image. From the bottom, the 
first three zeroes are ats = 0.5 + z'14.13, 0.5 + /20.02, 0.5 + i'25.01. Due to the relation to 
the Dirichlet eta function, the zeros also materialize at the same values of s, but on the 

q = 1/2 line. 

Remarkably, the blue streaks seem to be roughly parabolic, but are interrupted by nearly 
straight "erasures". The pattern is reminiscent of graphs of the strong-field Stark effect 
(need ref). In the Stark effect, eigenvalues are given by the characteristic values of the 
Mathieu functions. These cross over one another in a curious fashion; see for example, 

figure 20. 1 in Abramowitz & Stegun. 
The precise form of the parabolic blue streaks are due to a term of the form q~ s , as given 
by equation ll2.3l These may be subtracted, as illustrated in the next image. 



Figure .2. Hurwitz zeta, with leading term subtracted 




This image shows the magnitude of 

(,(s,q + l) = ^(s,q)-q- s 

for the range of < q < 1 along the horizontal axis, and < x < 50 along the vertical axis, 
just as in the previous image; the same coloration scheme is used as well. This essentially 
corresponds to the previous picture, with the leading divergence subtracted; for large 

positive x, the magnitude of the periodic zeta and the Hurwitz zeta differ by an 
exponentially small amount. A careful comparison of this picture to the previous now 
shows that the "erasures" or "streaks" in the previous image correspond exactly to the 
dark parts of this image. Again, a criss-cross lattice pattern can be seen. The dominant 
streaks in this picture are presumably due to the (1 + q)~ s term. 
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Figure . 3 . Extended view of Hurwitz zeta 




This figure shows a view of the magnitude of the Hurwitz zeta over an extended 

range, together with a rescaling of coordinates in an attempt to straighten out the 
striations. The value of q ranges from to 2\/2 « 2.8 along the horizontal axis, although 
it is scaled as g 3 / 2 along this axis. This rescaling essentially turns the parabolic striations 
into straight, radial lines emanating from the origin at the lower left. This image also 
re-scales the x coordinate, in an attempt to turn the radial striations into horizontals. 
The value of q = 1 is a vertical line running exactly down the middle of the image; the 
non-trivial Riemann zeroes are visibly stacked up on this line. The value of x runs from 
to 100 on this vertical line. On other verticals, x has been rescaled by q 1 ^ 2 , so that x runs 
from to 141 on the right hand side of the image, while being essentially held at zero 

along the left edge of the image. 
The increasing complexity of the chain-link pattern is apparent as one moves up the 
imaginary axis. Equally apparent is that the maximum range of excursion of the 
magnitude changes only slowly over the range of the image: as the pattern increases in 
complexity, it does not wash out, but has a distinct binary on/off character. 
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Figure .4. Poly logarithm phase plot 




This image shows a phase plot for the polylogarithm, on the complex z-plane, for 
s = 0.5 + /80. The color scheme is such that black represents the phase argLi s (z) = — 71, 
red represents the value of +7C, and green a phase of zero. Points where the phase wraps 
around the full color range are zeros. Thus, each of the "fingers" in this image terminates 
on a zero. Most of these zeros are close to, but usually not quite on the circle \z\ — 1. The 
lone zero near the center of the image corresponds to the zero at z = 0. The slight cusp at 

the right side of the image corresponds to z = 1 . The outer boundary of the image, 
adjoining to the solid green area, represents the limit of convergence for the combination 

of the Borwein algorithm and the duplication formula; one can enlarge the region of 
convergence slightly, but not much, but using polynomial approximations of higher order. 

The boundary is easily crossed by applying the inversion formula, as the later images 
show. The image is offset somewhat, rather than being centered on z = 0, because most of 
the convergent region is to the left of the imaginary axis. 
Note the other lone zero floating along inside of the |z| < 1 disk. If this picture were 
animated as T increased from in the positive direction, then the zero fingers can be see 

be seen as marching across in a counter-clockwise fashion, starting at z = land 
proceeding around. Most zeros pass to another sheet upon crossing z = 1, after making a 
full revolution; others spiral around a seconS time on this sheet, such as the lone zero up 
top. An animated movie of this image, showing the spiraling, can be seen at 
http ://w w w. linas . org/art-gallery/polylog/poly log .html 



Figure .5. Polylog whole plane 




This image illustrates the results of combining the Borwein algorithm together with the 
inversion and duplication formulas. The image shows the complex z-plane, is centered on 
z = 0, and runs from -3.5 to +3.5, left to right, top to bottom. The value of s is fixed at 
s = 0.5 + 15/ for the whole picture. The color scheme is the same as for the previous 

picture. 

Besides the zero at z = at the center of the picture, another zero is visible to the left, and 
slightly down. This zero is an "ex-Riemann zero", in that, if instead one had created the 

image for s = 0.5 + z'14. 1347 . . ., then this point would have been located exactly at 
Z = — 1. As it happens, this point has now rotated to slightly below z = — 1. The zeroes 
above and to the right of the origin are Riemann zeroes to be: as x increases, each will in 
turn rotate counterclockwise to assume the position at z = — I. 
The concentric shells are centered on the point z = 1 . The branch cut can be seen 
extending to the right from from z = 1, on the positive real axis. 
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Figure .6. Polylogarithm off the critical axis 




This figure shows the polylogarithm in the complex z-plane, at fixed 
s = a + iz = 1.2 + z'14. It is superficially similar to the previous image, with two notable 
differences. Visible just under the cut, on the right, is a zero. If a were smaller, this zero 
would move to the left; at a = 0.5, it would be very near to the branch point at z = 1. As x 
gets larger, this point moves up, crossing the branch cut and moving onto the next sheet. 

Of course, if instead one had (7 + z'x = 0.5 + /14. 13 . . ., this zero would be located 
precisely at z = 1 . Thus, this figure illustrates what a Riemann zero looks like when it is 

not on the critical strip. 
As x increases, the remaining zeroes circle around in a widening spiral: thus, the zero at 
the far left of the image is not near z — — 1, but is to the left of there (and thus is not going 
to be a zero of the Dirichlet eta function). If instead one had o with a value of less than 
0.5, the zeroes would spiral around in an ever-tightening spiral, and would never hit the 
branch cut or cross over. The Riemann hypothesis is essentially that these zeroes are 
impacting exactly at the branch point; a violation of RH would be a pair of zeroes that 
failed to hit the branch point, instead passing to the left and right of the branch point. This 
figure suggests that RH holds: the zeroes are clearly lined up in single file; it is hard to 
imagine how this single file might break ranks into a pair of zeros that simultaneously 

passed near trM branch point. 



Figure .7. Upper sheet of the polylogarithm 




This figure shows the g\Li s (z) = Li s (z) — Ai sheet of the polylogarithm on the complex 
z-plane, obtained by circling the branch-point at z = 1 in the right-handed direction. Two 
zeroes are visible: one at z — 0, and the other an "ex-Riemann zero", located not far from 



z = 1, a bit above the real axis. For this image, s is held constant, at s = 0.5 + z'15; the 
Riemann zeta function has its first non-trivial zero at s = 0.5 + z'14. 13 . . .. At that value, 
the above zero would have been located precisely at z = 1 , on the branch point of the 
polylogarithm. As x increases, the Riemann zeroes pass precisely through this branch 
point, leaving the first sheet and moving to this one. 
Indeed, the Riemann hypothesis can be taken as being equivalent to the statement that the 
Riemann zeroes always pass through the branch point in moving from one sheet to 
another. Taking s = a+ z'x, if a is held constant at a < 1/2 while x is increased, the 
polylogarithm zeroes never pass through the branch cut, and stay on the main sheet; 
whereas for a > 1/2, they do pass through the branch cut, but not at the branch point. 
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Figure . 8 . Lower sheet of the polylogarithm 




This figure shows the g j~ Li s (z) — Li s (z) — A_ i sheet of the polylogarithm on the complex 
z-plane, obtained by going around the branch-point at z = 1 in the left-handed (clockwise) 
direction. The concentric curves are centered onz= 1, and whose form is essentially that 
of (1 — q) l ~ s - The vertex of the black triangle approaches z — 0. By visual inspection, it is 
clear how to glue this sheet to the principal sheet shown in figure[3] After this gluing, a 
cut remains between the points z = and z = 1 . This cut may be glued to the sheet that 
results by winding around the branch at z = in a clockwise manner. The result is shown 

in the next figure. 
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Figure .9. More poly logarithm sheets 




Top row: shows the result of winding around z = 1, one, two and three times in a 
left-handed fashion. The cuts to the right of z =1 join together smoothly from image to 
image. Using the monodromy group notation, these are the g~[ , the g^ l g^ g{ , and the 

gi l g VrVo l 8i l sheets - 
Second row: Result of winding left-handed around z = 1, followed by winding around 

z — once. The cut < z < 1 joins smoothly to the images in the row above. The three 

images appear to be visually indistinguishable; but in fact they are not the same. Using 

the monodromy group notation, these are the gog\ l , the gog^gn l 8i 1 > an d the 

gogi l g l gi l g l gi l sheets. 
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